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Abstract. Monte Carlo techniques have been used to evaluate the statistical and systematic 
uncertainties in the helium abundances derived from extragalactic H II regions. The helium 
abundance is sensitive to several physical parameters associated with the H II region. In 
this work, we introduce Markov Chain Monte Carlo (MCMC) methods to efficiently explore 
the parameter space and determine the helium abundance, the physical parameters, and the 
uncertainties derived from observations of metal poor nebulae. Experiments with synthetic 
data show that the MCMC method is superior to previous implementations (based on flux 
perturbation) in that it is not affected by biases due to non-physical parameter space. The 
MCMC analysis allows a detailed exploration of degeneracies, and, in particular, a false 
minimum that occurs at large values of optical depth in the He I emission lines. We demon- 
strate that introducing the electron temperature derived from the [O III] emission lines as a 
prior, in a very conservative manner, produces negligible bias and effectively eliminates the 
false minima occurring at large optical depth. We perform a frequentist analysis on data 
from several "high quality" systems. Likelihood plots illustrate degeneracies, asymmetries, 
and limits of the determination. In agreement with previous work, we find relatively large 
systematic errors, limiting the precision of the primordial helium abundance for currently 
available spectra. 
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1 Introduction 

Standard big bang nucleosynthesis (SBBN) using the baryon density determined by WMAP 
[1, 2] predicts the initial abundances of D, '^He, ^He, and '^Li [3-10], allowing one to probe 
the early universe at redshifts of order lO^'' [11-13]. Therefore, the observed abundances 
provide a valuable check on the theory of SBBN, its concordance with the measurements of 
the microwave background radiation, and the content and interactions of the universe during 
the period of BBN [14-16]. To test these predictions, the observed abundances must be 
determined with relatively high precision. Because of the logarithmic relationship between 
the baryon to photon ratio, t], and the primordial helium abundance, Yp, the uncertainty 
of Yp must be < 1% to meaningfully test the theory. The 7- year WMAP value for rj is 
(6.19 ±0.15) X 10"i°, Komatsu et al. [2]. For comparison, the SBBN calculation of Cyburt 
et al. [10], assuming the WMAP r/ and a neutron mean life of 885.7 it 0.8 s [17], yields 
Yp = 0.2487 ± 0.0002, a relative uncertainty of only 0.08%. 

The determination of Yp is facilitated through low metallicity H II regions in dwarf 
galaxies. By fitting the helium abundance versus metallicity, one can extrapolate back to 
very low metallicity, corresponding to the primordial helium abundance [18]. The oxygen 
to hydrogen ratio, O/H, commonly serves as a proxy for metallicity. Though this area of 
research has benefited from three decades of development, the determinations of Yp have 
suffered from significant differences between the results. The difficulties in calculating an 
accurate and precise measure of the primordial helium abundance are well established [19- 
21]. Here, we introduce a new method based on Markov Chain Monte Carlo (MCMC) 
techniques. 

Observations of the helium to hydrogen emission line ratio from extragalactic H II 
regions provide a measure of the helium to hydrogen ratio, y^ = -^^j-jjy- Correspondingly, 
the statistical measurement errors in the helium and hydrogen emission line fluxes contribute 
to the uncertainty on y"*". Unfortunately, this calculation of y^ and its uncertainty are 
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complicated by a myriad of systematic effects. Interstellar reddening, underlying stellar 
absorption, radiative transfer, and collisional corrections alter the observed flux, complicating 
the measurement of y"*", and amplifying the uncertainty. The photons are scattered by dust on 
their journey (interstellar reddening). The stellar continuum juxtaposes absorption features 
under nebular emission lines (helium and hydrogen underlying absorption). The H II region 
itself absorbs and re-emits photons (radiative transfer); both recombination and collisional 
excitation contribute to the emission (collisional corrections for helium and hydrogen). None 
of these processes can be directly measured and, therefore, cannot be removed independent of 
the observed emission lines and theoretical models. As a result, the uncertainty on y"*" must 
reflect the presence of, and lack of certainty regarding, these systematic effects. Determining 
y"*" in conjunction with the robust estimation of the model parameters used in correcting for 
the listed systematic effects requires "high quality" spectra. This desired confidence weighs 
against the need for larger sample sizes to decrease the uncertainty on Yp (and dY/dZ). 

The importance of Monte Carlo techniques was demonstrated in a "self-consistent" 
analysis method, stemming from the work of Izotov, Thuan, & Lipovetsky [22] and Peimbert, 
Peimbert, & Ruiz [23], for determining the nebular helium abundance based upon six helium 
and three hydrogen lines [19, 20]. In preceding work [24], hereafter AOS, we updated and 
extended the physical model and integrated the helium and hydrogen calculations with the 
goals of improving accuracy and removing assumptions. The focus of the current paper is 
the exploration of a new technique based on a Markov Chain Monte Carlo (MCMC) analysis 
and departs from the "self-consistent" method. Rather than fitting the parametric inputs 
to a helium abundance, the frequentist approach developed here builds a global likelihood 
function for all parameters including the helium abundance. 

As we will demonstrate, the MCMC method is statistically superior to previous efforts, 
it is more direct and transparent, and it maintains efficiency. Of primary benefit, the results 
are more rigorous: the solution remains unbiased by the procedure, the uncertainty captures 
the confidence of the model and measurements, visualization of the parameter space topol- 
ogy reveals the reliability of the determination, and spectra failing to resolve their physical 
environments are identified. 

Section 2 discusses the determination of parameter uncertainties and details the differ- 
ences in our approach between AOS and this work. The utility of Markov Chains and the 
computational implementation are described in §3. In §4, we describe tests using synthetic 
data to illustrate the method and its utility, with particular emphasis on secondary min- 
ima and the incorporation of a temperature prior. MCMC is implemented in analyzing the 
dataset used in AOS in §5, and, in §6, Yp is determined. Finally, §7 offers a discussion of the 
exploration and results as well as the next steps in better determining the primordial helium 
abundance. 

2 Parameter simulation 

The purpose of the Monte Carlo is to generate a set of parameters with their values. 
Improving upon the previous efforts, a new technique, MCMC (see §3), enables simultane- 
ous Monte Carlo over all model parameters including the He abundance. A full frequentist 
sampling of the parameter space involving eight input parameters would be computationally 
prohibitive. However, using MCMC allows one to judicially sample regions of parameter 
space with relatively low x^- In this way, we can construct the likelihood function and deter- 
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mine the best fit point in the multidimensional parameter space along with their associated 
uncertainties. Indeed, the uncertainties are the primary focus of this work. 

The Monte Carlo approach of ref. [20] and AOS took each set of measured fluxes and 
built a Gaussian distributed dataset of fluxes based upon their measurement uncertainty. For 
each of a 1000 such datasets, a best-fit solution was found for the helium abundance as well as 
the physical input parameters using the "self-consistent" method which determines the set of 
input parameters with a based on the derived helium abundance from each of six helium 
emission lines. The final result was computed from the average and standard deviation of the 
set of solutions. Using the fiuctuation of the minimum is, however, not a direct measure of the 
X^'s parameter dependence. Furthermore, it is also not as robust as desired. Each solution 
was restricted to physically meaningful parameter space (e.g., positive densities), potentially 
biasing the solution. Additionally, as was manifested in AOS and will be discussed further 
in §4, functions lacking a well constrained temperature and density can produce unlikely 
high density and low temperature solutions that greatly skew the results. Ultimately, these 
considerations, tempered by the required computational efficiency, motivate this work. 

The x^ function defined here, and used for parameter fitting, is modified from that 
used in previous work. Rather than defining y"*" implicitly, as the average of six individual 
line abundances, and minimizing the deviation between the lines, y"*" is demoted to an input 
parameter, no different than the others (e.g., temperature and density). Instead, here, we 
use all of the input parameters (described below) and calculate synthetic fiuxes which are 
then compared to observed flux, weighted by the observed uncertainty, allowing for a more 
standard definition of x^, 

( ^(A) F{X) .2 

meas (n ^\ 



where the He flux at each wavelength A relative to the fiux in H(3 is given by 
F(A) _ + L{X) w{Hi3) „ ... 1 + 7j(A) 



fr{X) 'V 10-/W^(^/^). (2.2) 



W{X) 

The x^ in eq. 2.1 runs over He and H lines, and the ratio of H fiuxes is defined analogously, 

^ ^(A) W{HI3) l + 7;(^) .r.-f(X)C{Hl5) ON 

For the above flux equations, six measured helium emission line fluxes (A3889, 4026, 4471, 
5876, 6678, and 7065) and three hydrogen emission line fluxes {Ha, H-f, H6), each relative 
to H/3 ( _p(gff ) ) ; along with their equivalent widths {W{X)) are used. The predicted model 
fluxes are calculated from an input value of y"*" and emissivity ratio of H/3 to the helium or 
hydrogen line, ^^x) ' '^ith corrections made for reddening (C(H/3)), underlying absorption 
{a,H <fe a/fe); collisional enhancement, and radiative transfer. The optical depth function, 
fr, and collisional to recombination emission ratio, ^, are both temperature (T) and density 
(ug) dependent (the emissivities are also temperature dependent). Additionally, the hydrogen 
collisional emission depends on the neutral to ionized hydrogen ratio (^). Therefore, there are 
a total of eight model parameters (y+, ng, a^e, 7", T, C(H/3), a//, ^). The physical model itself, 
the equations relating the abundance and correction parameters to the flux, is unchanged from 
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AOS, and incorporates all of the effects explored therein (updated emissivities, wavelength 
dependent underlying absorption, and neutral hydrogen collisional emission). 

We note here that the equivalent width of an emission line is not independent of the 
flux. Each flux is related to its equivalent width through the flux of the continuum at that 
wavelength, h(A), as 

F{\) = W{\)h{X). (2.4) 



As a result, if the model parameters generate a lower flux, the equivalent width should be 

lowered correspondingly. The continuum flux is constrained to remain constant such that 
fe(A) 



the ratio -rrwL: is determined from the measured flux ratio and equivalent widths. 



hiX) F(A) W{Hl3)^eas 

(2.5) 



hiHP) F{H(3)„ W{\)meas ' 

Eq. 2.2 can therefore be rewritten to remove W{\) entirely and solve for a consistent emission 
line ratio to H(/3). This yields a simplified He flux ratio, 

F{\) + E{X) l + g(A) f^,)c(H0) W{H/3) + aH{H/3) UHeW h{X) 

F{HI3) y E{Hp)^^^ 'l + c^Hf3) W{Hp) W{Hp)h{Hpy 

(2.6) 

and H flux ratio, 

F{\) ^ E{\) 1 + %{X) f^,^c(H0) WiHl3) + aH{H/3) _ aniX) hjX) 
F{Hp) ^(F/3) l + g(F/3) W{H^) W{HI3)h{Hpy ^'^ 

Because it is this ratio of fluxes that is calculated, W{Hj3) cannot be removed from the 
equation. The treatment of this measured quantity (equivalent to the H/3 flux measurement) 
will be discussed in the next section (§3). 

To reduce bias in the solutions of AOS, the emissivities of [25, PFM] were refit with a 
stronger density dependence to protect the minimizations from running away to unphysical 
regions of the parameter space. As mentioned above, the method of AOS was susceptible 
to the skewing of the final average. Because at large densities the density dependence of 
the PFM emissivities is negligible, extreme outliers with densities increasing nearly without 
bound occurred, rendering the final result meaningless. Refitting the emissivities to the form 
of Benjamin, Skillman, & Smits [26] rectified this behavior, while preserving the accuracy 
of the PFM equations in the region of interest. However, the approach of this work, direct 
Monte Carlo over the parameters, is not susceptible to the previously exhibited biasing. A 
poor constraint on the density will be manifested in the function. This will affect the 
uncertainty in a transparent way, reflecting the quality of the determination. The solution, 
the minimum x^5 will be unaffected, and any degeneracy in the solution, i.e., multiple local 
minima, will be apparent. As a result, the unmodified PFM emissivities [25] are used in this 
work. 



3 Markov Chain Monte Carlo 



The Markov Chain Monte Carlo (MCMC) method is an algorithmic procedure for sampling 
from a statistical distribution [27, 28]. Therefore, the sequence of points in the parameter 
space reconstructs the target distribution, here the likelihood distribution {C). The real value 
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of MCMC is the judicious choice of samphng such that the density of samples is greatest 
around the best-fit point and increasingly sparse in increasingly unlikely regions in parameter 
space. As a result, computational time is not wasted in exploring poor model solutions, while 
full exploration of complex topology is still facilitated. In brief, the algorithm is efficient, with 
computational time scaling roughly linearly with the number of parameters. This is in sharp 
contrast to the exponential growth in time cost encountered with gridding. For example, if 
only 100 points in each parameter were taken, then for eight parameters, a total of (10^) 
evaluations would be needed. Typically, the results of this work are based on 10^ — 10^ points. 

CosmoMC^ is a freely available Fortran package for MCMC, designed originally for 
WMAP parameter extraction, but easily modified for generic sampling. It allows for a 
variety of sampling algorithms, but this work makes use of the classic, and widely applicable, 
Metropolis-Hastings algorithm [28, 29]. We must also choose a proposal function which is a 
distribution function determining the sizes of the steps to take in each parameter. The most 
common choice, and the one used here, is an n-dimensional Gaussian distribution, where n 
is the number of input parameters. For a symmetric proposal function, with points in the 

Markov chain: 

• Select an initial parameter vector, xq 

• For i = to A - 1: 

1. Randomly choose Ax from a proposal function 

2. Xj+i = Xi + Ax 



4. Accept Xj+i with probability r (if r > 1, r = 1) 

5. Otherwise Xj+i = Xi 

In effect, a Gaussian distributed random walk is used; whereby, points that are of higher 
likelihood are always accepted, and points of lower likelihood are accepted proportional to 
the ratio of the proposed to current likelihood. For this work, the likelihood is simply, 

C = exp(-xV2), (3.1) 

with given by eq. 2.1. After a 'burn- in' period, the points will generate a Markov Chain 
converging to the target distribution. The most important tuning comes in picking the 
proposal function step size, the Gaussian width. If this is too large, then the acceptance 
fraction will be very small, and the points will be clumped with large separations. If the 
width is too small, the exploration of the full distribution will be very slow, and failure to 
explore all minima becomes of increasing concern. Roughly, a proposal width similar to the 
target distribution width, for each parameter, is the most efficient. Furthermore, proposal 
widths leading to an acceptance fraction of ~25% have been shown to be optimally efficient 
[30, 31]. 

One complication in this analysis comes from the emission line equivalent width of H/3, 
W{Hf3). It is used in making the corrections for underlying absorption, and is, therefore, 
needed to calculate the predicted flux, even though it is itself a measured quantity. To 
treat this, W{HP) is encoded as a nuisance parameter. For each MCMC point, a new 
Gaussian distributed equivalent width is chosen based upon the measured value and error. 

^http: //cosmologist . info/ cosmomc/ 
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A corresponding term is then also added into the x^j 



2 _ {W{HP)-W{HPU,as? 

" '^{m^ ■ ^"^-^^ 

FinaUy, the use of priors is straightforward in this MCMC analysis. The restriction 
of the parameters to physically meaningful values, such as densities greater than zero, just 
means that proposal points in unphysical regions are automatically rejected. Additional mea- 
surement information is incorporated directly into the likelihood in the form of an additional 
contribution. For example, a temperature characterizing the H II region and derived from 
other emission lines could be incorporated as, 

C = exp(-xV2) exp(-(T - Tmeas? I'ia^)- (3.3) 

Consequently, points with temperatures diverging significantly from T^eas would be disfa- 
vored by their low likelihood, and the uncertainty on the best-fit temperature, and potentially 
all other parameters, would be reduced. 



4 Synthetic exploration 

4.1 Improved analysis with MCMC 

Synthetic data are useful to demonstrate the utility and simplicity of the MCMC analysis. 
Synthetic helium and hydrogen fluxes were generated using characteristic parameter values 
of y+ = 0.08, Ue = 100 cm-3, am = 1-0 A, r = 0.2, T = 18,000 K, C(H/3) = 0.1, a.H = 1-0 A, 
and ^ = 1.0 X 10~^. Flux errors were 1% for the hydrogen and 2% for the helium lines with 
EW(H/3) = 250 A. Given the synthetic fluxes generated from these input parameters, we 
can apply the MCMC routine to determine the likelihood function across our 8-dimensional 
parameter space. Figure 1 shows the ID marginalized likelihood for the helium abundance. 
Each point shown represents a minimization of the function over variations of the other 
seven input parameters at a given (binned) value of y"*". As one can see, the distribution 
nicely reproduces the input value of y"*" with the 68% CL range indicated by the dashed 
lines. Figures 2 and 3 show the other seven parameters. The 68% CL ranges for all of the 
parameters are collected in table 1. 

The likelihood plots are revealing in several important ways. First, the similarities 
between the abundance and temperature are not a coincidence (figures 1 and 2d). The 
temperature is the most important parameter in determining y+ through the emissivity 
equations (see, e.g., AOS). Also, because of a degeneracy in the determination, primarily 
between the temperature and the density (cf. AOS), the temperature is not well constrained 
and admits values 7,000 K below the input value at the 95% level. As shown in figure 4, 
there is a strong negative correlation between these parameters which allows a large variation 
in their values without rapidly raising . The 95% levels, defined by the x^ = 4 contour, 
allow densities reaching 1250 cm~^ and temperatures ranging between 11,000 K and 22,000 
K, and result in the abundances spanning 0.072 to 0.089. 

Second, the neutral to ionized hydrogen ratio (figure 3c) is unconstrained at large val- 
ues, exhibiting asymptotic behavior at x^ = 4. The collisional relative to recombination 
contribution to the hydrogen emission depends exponentially on the temperature (see AOS), 
and as a result, the low temperatures arising out of the temperature-density trade-off permit 
nearly any value of the neutral hydrogen fraction, including unphysical ones. In the previous 
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Figure 1. x versus abundance for synthetic data with model parameters, y 
aue = 1-0 A, T = 0.2, T = 18,000 K, C(77/3) = 0.1 
confidence level is marked by the dashed lines (x^in 
the input value. 



0.08, rie = 100 
an = 1.0 A, and ^ = 1-0 x 10"^ The 68% 
= 0.0 for synthetic data). The arrow denotes 



approach, this behavior lead to some of the flux perturbed solutions with very large neutral 
hydrogen fractions. This greatly biases the final averaged result, ^ = 27.24 x 10~^, away from 
the input value of 1.0 xlO~^ (table 1). This biasing is exacerbated by the restriction of solu- 
tions to physical values, prohibiting negative solutions within the dataset. Such restrictions 
are physically meaningful and necessary but, because of the statistical approach, impact the 
final result detrimentally. In less pronounced fashion than for the neutral hydrogen fraction, 
any parameter with an unperturbed best-fit point near zero suffers from this effect in the 
analysis of AOS. The optical depth, r, with a generating value of 0.2, exhibits this biasing 
clearly, returning 0.42. For the singular solution of the new MCMC analysis, the restriction 
to physical values is correct and carries no possibility of unintended bias. 

The results based on the method of AOS are also shown in table 1 and can be directly 
compared with the current results based on MCMC. The AOS result reproduces the input 
helium abundance fairly well (y"*" is 1.4% high which is less than one third of the 4.3% 
uncertainty, see table 1). In contrast, using MCMC, the final result for the abundance is that 
of the true minimum, and any asymmetry in the likelihood impacts only the uncertainty. 
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Figure 2. Similar to figure 1 with plots of versus density, helium absorption, optical depth, and 
temperature. 



The same is true for the other parameters as welh 
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AOS 


MCMC 


He+/H+ 

He 

ABS(He I) 


0.08 
100.0 
1.0 


0.08108 ± 0.00341 
147.8 ± 282.8 
1.05 ± 0.09 


07999 +^-^^''^^^ 
u.u^yyy -0.00292 

100 1 +^^^-^ 
1 00 +°-°^ 

^■^^ -0.06 


T 


0.2 


0.42 ± 0.41 


20 


elm 


18,000 
0.1 


17,440 ± 2308 
0.08 ± 0.03 


1 7 qqq +2239 
± ( -2711 

in 

^■-■^^ -0.03 


ABS(H I) 


1.0 


1.35 ± 0.84 


1 00 

J^-^U -0.64 


e X 10^ 


1.0 


27.24 ± 187.27 


98 +59-77 
-0.98 



Table 1. Comparing Gaussian Distributed Fluxes (from AOS) and MCMC Analyses with Synthetic 
Data 



-8- 



C(Hp) 
0.06 0.10 




500 1000 1500 2000 2500 
Ux10') 

Figure 3. Similar to figure 1 witlr plots of ^ versus reddening, hydrogen absorption, and neutral 
hydrogen fraetion. 



4.2 Degeneracies at large values of optical depth 

An important feature emerges in synthetic testing at high optical depth. Using the same 
generating values as in §4, except for r = 2.0 instead of 0.2, results in a pronounced second 
minimum in the helium abundance (see figure 5). For objects with significant radiative 
transfer and relatively large errors on He A3889, the parameter space contains the flexibility 
to decrease r abruptly from the best-fit value, at a correspondingly low temperature, low 
abundance, and high density, to decrease x^- Consequently, a second minimum is found in 
the MCMC analysis. 

The reason for this second minimum can be understood by examining figure 5 in ref. [19]. 
The He I A7065 emission line has a strong dependence on density (it is increased via collisional 
excitations from the meta-stable triplet 2s level). However, at significant values of optical 
depth, the He I A7065 emission line is also increased by absorptions out of the meta-stable 
triplet 2s level which cascade down through this emission line. Because the He I A3889, is also 
increased by via collisional excitations from the meta-stable triplet 2s level, but decreased by 
absorptions from that level, the two lines can trade against each other, combined with the 
temperature-density degeneracy, to form the second, non-physical, minimum. 

In figure 5, the second local minimum is near y"*" = 0.072, and correspondingly, T = 
10,000 K, Ug = 1500 cm~^, and r = 0.75. Though it does not impact the 68% CL, it lies well 
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Figure 4. Contour plot of versus density and temperature for tlie same synthetic model used in 
figures 1-3. Contours are Ax^ = {1, 2.3, 4, 6, 9}. The degeneracy between the parameters limits the 
precision of their determination and, therefore, of the helium abundance. 

within the 95% CL limit {x^=l.6). Results fr om the previous method (AOS) were in fact 
susceptible to being skewed by secondary minima. Upon perturbation of the fluxes, this lower 
minimum could become the best-fit point for some of the datasets. Then, upon averaging, 
these outlier solutions bias the final result. As already discussed, the MCMC solution is 
unaffected by asymmetries about the minimum; however, broad, deep secondary minimum 
will greatly, and artificially, increase the calculated uncertainty. Furthermore, for real data, 
the presence of a second local minimum, with a similar x^, raises doubts regarding which 
minimum represents the physical environment. The object's reliability is clearly decreased, 
but, in the following section (§4.3), the use of an additional observation of the temperature 
via the [O III] emission lines to resolve this ambiguity is discussed. 

This work uses the radiative transfer calculations of Benjamin, Skillman, &: Smits [32] 
[first used in 20] . These fits assume a spherical nebula with uniform density and temperature 
and no systematic velocity gradients. However, H II regions may be expanding - the original 
work of Robbins [33] models the regions as expanding with a constant velocity gradient. 
Similarly, the effects of turbulence may be more complicated than as incorporated. Any such 
deviations from the model will be more pronounced at larger optical depth. Furthermore, 
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Figure 5. versus abundance for synthetic data with model parameters, ~ 0.08, rig ~ 100 cm ^, 
aHe = 1-0 A, T = 18,000 K, C(iJ^) = 0.1, an = 1.0 A, and ^ = 1.0 x lO"''. The only difference 
between the two sets is the value of the optical depth; r ~ 0.2 for the lighter, open squares and r = 2.0 
for the darker, solid circles. For the latter, the 68% confidence level is marked by the dashed lines, 
and the input value is denoted by the arrow. The larger optical depth allows a secondary minimum 
to develop at low abundance. Such behavior highlights a deficiency in the model and mitigates the 
reliability of galaxies exhibiting large optical depth. 



He A3889 is the primary line in determining r; yet, it is blended with H8, decreasing its 
reliability. The combination of these concerns, in addition to the emergence of the secondary 
minimum, motivate a preference for objects with low optical depth. 

4.3 Using T(0 III) as a conservative prior 

Adding an independent measurement of the temperature in the H II region has the potential 
to weaken the temperature-density degeneracy. The [O III] AA4363, 4959, and 5007 emission 
lines provide such and independent measurement, T(0 III), associated with the doubly ion- 
ized oxygen of the H II region. However, because the temperature in the H II region is not 
expected to be perfectly uniform, the exponential sensitivity to temperature of the [O HI] 
emission lines biases the derived temperature to higher than average values. Thus, T(0 HI) 
may not well characterize the average electron temperature over the entire He II emission 
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zone. The magnitude, and thus importance, of this temperature difference is a matter of de- 
bate. Recent work by Peimbert, Peimbert, & Luridiana [34] found that T(He II) was always 
less than T(0 III), ranging between 6-10% less (solved) and 3-11% less (photo-ionization 
models). However, from comparisons of temperatures derived from the Balmer jump with 
temperatures derived from the [O III] emission lines, Guseva et al. [35, 36] find no evidence 
for this temperature offset from studies of low metallicity H II regions. 

To use T(0 III) directly as a measurement of the average temperature would break the 
temperature-density degeneracy (with resulting very small uncertainties), but risks biasing 
the derived values of helium abundance to artificially high values in the presence of non- 
uniform temperatures. In fact, T(0 III) can be viewed as a limit to the upper bound 
on the average temperature in the H II region, and as a consistency check to eliminate 
physically improbable low temperature solutions. The goal is to use T(0 III) to constrain 
T(He II) to physically meaningful values, not to bias it. Note that the expected differences 
between T(0 III) and T(He II) are much larger than the calculated error on T(0 III). 
Therefore, in incorporating T(0 III) as a prior, this work takes (t{Toiii) = 0.2 Tom, a 
very weak constraint. This will keep bias in the solution negligible, while employing the 
available additional temperature diagnostic to better isolate the physically relevant parameter 
space. In the case of likelihoods exhibiting a second minimum with a very low corresponding 
temperature (as demonstrated in §4.2), the T(0 III) prior rigorously and definitively rules 
out the unphysical region. As a result, the parameter errors are then well defined, and the 
solution is unambiguous. 

Figure 6 illustrates the effect of using the conservative T(0 III) prior on the synthetic 
object's likelihood. The low temperature local minimum is completely removed by using the 
T(0 III) measurement as a prior, even with the very gentle implementation chosen. As ex- 
pected, due to the temperature-density degeneracy, the increased constraint on temperature 
translates into a marked improvement in the derived density, as is evidenced in figure 7. The 
only concern is the possibility of bias in the best-fit solution. For this example, T(0 III) 
was input as 19,000 K, 6% higher than T(IIe II). Table 2 compares the minimum solution 
with and without the use of the prior. The increased temperature leads to y"^ = 0.0802 (at 
= 0.05), but this is only 0.3% high compared with an uncertainty of nearly 5%. This 
uncertainty (on y"*") itself benefits from the prior, especially on the lower side, as would be 
expected due to Tqiii > ThcII- In effect, the inclusion of the [O III] temperature prior, 
via our conservative implementation, produces a negligible bias, reforms parameter behavior, 
and strengthens the determination of the solution. 

5 Application of MCMC to real observations 

To illustrate the effect of using the MCMC method, we will calculate the resulting value 
of the helium abundance y"*", along with the other input parameters and their associated 
uncertainties, for the seven galaxies from Izotov & Thuan [37] analyzed in ref. [20] and AOS. 
Additionally, NGC 346 [23] and I Zw 18SE [38] are also included, as in AOS. The values of 
T(0 III) used are taken from ref. [20]. 

MCMC, with the T(0 III) prior, was used to determine the likelihood function for 
each of the nine objects considered in AOS. The best fit value of the input parameters and 
their uncertainties are found from the likelihood function as described in section §4. Mrk 193, 
SBS 1420+544, and SBS 0335-052E each exhibit a pronounced double minimum. As expected 
from the synthetic testing in §4.2, these three objects have the largest optical depths. 
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Figure 6. versus abundance for synthetic data with model parameters, — 0.08, rie — 100 cm ^, 
aHe = 1-0 A, r = 2.0, T = 18, 000 K, C{HP) = 0.1, an = 1.0 A, and ^ = 1.0 x lO"'^. The darker, solid 
circles are analyzed with T{0 III) — 19, OOOA' included as a prior, while the lighter, open squares do 
not include the prior. The 68% confidence level is marked by the dashed lines, and the input value is 
denoted by the arrow. Note that the minimum is not noticeably impacted by the prior. 
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Table 2. Comparison of the effect of a T(0 111) prior with Synthetic Data 
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Figure 7. Contour plot of versus density and temperature for ttie same synthetic model used in 
figure 6. Contours are Ax^ = {f , 2.3, 4, 6, 9}. The dramatic curtailing of the density and temperature 
variance is due solely to the inclusion of the weak T(0 III) prior - with and without in the top and 
bottom panels, respectively. 
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As an example of the degeneracies for the large optical depth regime, we examine 
Mrk 193. Shown in figure 8, the T(0 III) prior clearly distinguishes the physical from 
the unphysical minimum. Before including the O III measurement, the two minima are 
troublingly similar in their values, yet dramatically different in their physical environment. 
The slightly higher minimum with ^ = 4.3, has a lower abundance, lower temperature, and 
higher density than the one with = 4.0. The parameter values of the higher minimum are 
y+ = 0.07271, Ue = 5675 cm~^, and T = 6,655 K. These values of density and temperature 
are physically unreasonable. The lower minimum has values, y^ = 0.08572, Ug = 2 cm~^, and 
T = 14,017 K, which are certainly quite reasonable. Although each minimum defines a 68% 
CL by itself, neither does so at the 95% CL. Therefore, even though the lower minimum is the 
physically relevant one, their near equivalence in and poorly constrained uncertainty values 
would undermine the reliability of this object. By including the [O III] temperature {Tqiii = 
15, 561 K), the unphysical minimum is completely removed, and the galaxy produces a well- 
defined solution: y+ = 0.08619, Ug = 1 cm~^, and T = 15,226 K. This new minimum is shifted 
0.5% higher, but this variance is insignificant in comparison with the overall uncertainty of 
5%. Unlike in the case of synthetic testing, the difference in the solutions with and without 
the T(0 III) prior is not a measure of the bias itself, but of an upper bound on the bias. 
Here, the solution excluding the [O III] measurement is not necessarily more accurate and, 
in fact, can be seen as less accurate due to its neglect of a relevant observation and larger 
uncertainty. 

The remaining six objects are much better behaved. Figure 9 exhibits the solution 
for one of these, SBS 1159+545. For these objects the shapes of the likelihood plots are 
similar to those found in the synthetic results of figures 1, 2, and 3, prior to the inclusion 
of the [O III] prior. For SBS 1159+545, the T(0 III) prior noticeably increases the helium 
abundance found in the solution. This occurs primarily due to the higher temperatures 
found when the prior is included (using the Tqiii value of 18,600 K for SBS 1159+545 
raises the solution temperature, T, from 15,900 K to 17,100 K). Reassuringly, this large 
temperature difference only raises the abundance by 1% (0.08326 to 0.08416), compared with 
an uncertainty of 5%. As was previously mentioned in discussing Mrk 193, the solution with 
T(0 III) is preferred because the added information strengthens the physical determination 
and reliability, outweighing the small unwanted bias. 

Figure 10 is included to illustrate the importance of using the T(0 III) prior conser- 
vatively. Though figure 8 clearly shows the dramatic benefit of incorporating the [O III] 
measurement, we chose a very conservative implementation. Had we used a more strict 
T(0 HI) prior, any possible deviation between T(He II) and T(0 III) would have been 
washed out. Interestingly, three objects show He temperatures greater than the T(0 III) 
temperature, though they are all within la of being equal. The remainder of the objects do 
show lower He temperatures, as might be expected, and it is significantly lower for several 
galaxies. Therefore, our weak prior still allows the solved temperature to refiect the physical 
environment defined by the helium emission. 

The results, best-fit solution and uncertainties, for the entire sample are presented 
in table 3. These are compared with the results from AOS and are found to be similar 
(differences of less than 1 a) for the all of the objects. The most significant difference in 
the helium abundance is for SBS 0335-052E, where the abundance likelihood has a very 
broad, shallow bottom and a weak second minimum. As a result, the previous solution 
for SBS 0335-052E was determined by an average of abundances that spanned this large 
range in y+, including the weak secondary minimum with a low abundance and temperature. 
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Figure 8. versus abundance for Mrk 193. The darker, solid circles are analyzed with the T(0 III) 
prior, while the lighter, open squares are without. The 68% confidence level is marked. The prominent 
double minima, though similar in x^, correspond to very different temperatures and, therefore, are 
readily resolved by including the T(0 III) prior. The shift in the minimum is negligible in comparison 
to its uncertainty. 



Therefore, the new solution, consistent with the T(0 III) prior, is much higher. In the three 
cases with the largest shifts in helium abundance, SBS 0335-052E, SBS 0940+544N, and 
NGC 2363A, the helium abundance increased due to an increase in the temperature, which 
results from the use of the T(0 III) prior. 
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Figure 9. versus abundance for SBS 1159+545. The darker, solid circles are analyzed with the 
T(0 III) prior, while the lighter, open squares are without. The 68% confidence level marked. Very 
similar to synthetic data (figure 1), this object is characteristic of well-behaved parameter determi- 
nations. The notable shift in the minimum is due to a large difference in T(0 III) and T(IIc II), but 
is still much less than the uncertainty. 
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Table 3. Comparison of Physical Conditions and Hc+/H+ Abundance Solutions 



6 Results of the MCMC analysis 



A comparison of the derived values of y"*" presented in table 3 is plotted in figure 11. An 
important result of this comparison, between the previous analysis of AOS and the MCMC 
method, is that the uncertainties are very similar. Indeed, the maximum uncertainty (the 
larger of the two asymmetric uncertainties) is larger but, in all cases, is much less than 1 a 
larger. On average, the maximum uncertainty increased by 17%. 

That the uncertainties increased at all might be unexpected. Naively, one would expect 
that the use of the conservative prior, which eliminates an unphysical minimum, would result 
in a decrease in the uncertainty. However, the comparison between the analysis in AOS and 
the MCMC method with prior reveals that the previous method inherently underestimated 
the errors. For example, for Mrk 193 (figure 8), the available parameter space for the Monte 
Carlo analysis was much broader for the AOS analysis, but the estimate of the uncertainty 
did not capture the true uncertainty. Primarily, this is because the variance was calculated 
symmetrically about the mean, not asymmetrically about the minimum, and because the 
method was not a rigorous approach to exploring the parameter variation. The MCMC 
method yields an improved statistical result and more accurately reflects the uncertainty of 
the helium abundance determination. 

The development of the new MCMC analysis technique is the primary goal of this work. 
To aid in demonstrating the effect, we also calculate the primordial helium abundance (mass 
fraction), Yp. A regression of Y, the helium mass fraction, versus 0/H, the oxygen to hy- 
drogen mass fraction, is used to extrapolate to the primordial value. ^ For direct comparison, 
the 0/H values are taken from AOS [which took the values from 20, for the same reason]. 
The relevant values are listed in table 4. 

NGC 346 is more chemically evolved, and I Zw 18SE has a low equivalent width of 
H/3, making it more sensitive to underlying absorption, and is a candidate for Galactic 
Na I absorption due to its radial velocity. As a result, neither object was included in the 
primary regression of AOS and will similarly be excluded here. This regression yields Yp = 
0.2609 ± 0.0117, with slope -67 ± 214 and = 1-7. Note that dY/dZ is very poorly 
determined and if we restrict the analysis to theoretically meaningful positive slopes, we find 

yp = 0.2573t™, (6.1) 

with a range in slopes from to 149 and a of 1.8. Equation 6.1 agrees with the WMAP 
value of 1^ = 0.2487 ± 0.0002 within la. We do note, however, that the two most recent 
measurements of the neutron mean life [39, 40] are significantly lower than the accepted world 
average [17]. A drop in the mean life of about 6 s would result in a lower BBN abundance by 
about 0.001. While, conservatively, it is premature to claim a discrepancy in helium (because 
the systematic uncertainties are so large), if real, it may hint towards new physics [41-43] or 
require an early astrophysical source for helium [44]. AOS determined Yp = 0.2561 it 0.0108 
(0.25611qq5o8 when the slope is restricted to be positive); that these results agree so closely 
is not surprising given that the dataset and physical model of the works are identical. 
As the 0/H domain is very limited, a mean evaluation is justified and gives, 

Yp = 0.2573 ± 0.0033. (6.2) 

Inclusion of I Zw 18SE and NGC 346, the lowest and highest metallicity objects, reduces the 
intercept and error to 1^ = 0.2549zb0.0066 with a slope of 45 ± 86. Much of the improvement 

^This work takes Z = 20[O/H) such that Y = 4y(i-20(o/g)) 
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Figure 11. Abundance comparison for tlie target objects as analyzed in AOS, using Gaussian 
distributed fluxes, and MCMC. The lighter, thinner lines are for Gaussian fluxes, with MCMC given 
in thicker, darker bars. The quoted primordial helium abundance, Yp, is based on a regression of 
the seven objects of AOS (i.e., I Zw 18SE and NGC 346 have not been used in the regression). This 
sample is also used to produce the mean helium abundance, < Y >. 



stems from the longer metallicity baseline but, therefore, also more assumptions. It has always 
been our view that one well measured object at high metallicity, which would fix dY/dZ, 
should not affect the uncertainties of measurements at low metallicities if their uncertainties 
are large. The size of our overall uncertainty also suffers from a low number of "high quality" 
measurements at low metallicity which should be close to primordial. Nonetheless, this is 
the state of current data. Including more "high quality" objects would clearly benefit the 
determination. This work has focused on introducing a new statistical method, saving a 
revised sample for later work. Indeed, we hope that this method will allow the unambiguous 
use of other low metallicity data which are found in the literature. 



7 Summary 

The primary result of this work is the demonstration of a statistically rigorous method for 
determining the uncertainty of the abundance and model parameters. The use of MCMC 
allows one to efficiently sample the parameter space so that the uncertainties can be calcu- 
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Table 4. Primordial Helium Regression Values 



lated directly from the change in as the parameters are varied from the best-fit solution. 
Computationally, MCMC is efficient and straightforward. Beyond the improvement in the 
approach itself, the constructed likelihood distributions for the parameters are instructive in 
evaluating the quality of the object and illuminating the sources of differences with previous 
analyses. 

Particularly illustrative of the benefits of the likelihood approach was the discovery of the 
increased degeneracy at large optical depth. With synthetic data, a second minimum emerges 
as the optical depth increases. As this false minimum becomes more significant, the reliability 
and quality of the object is undermined. In concordance with this predicted effect, several 
galaxies, each with large optical depth, exhibit prominent second minimums. However, a 
second benefit of the approach is the straightforward incorporation and interpretation of 
priors. This allows the electron temperature defined by the [O III] emission lines to be 
utilized to eliminate low temperature, unphysical secondary minimums. Therefore, after 
taking the [O III] measurement into account, large optical depth objects are well behaved. 
It also worthy of note that the prior is used very conservatively. This ensures that the solved 
value of the temperature primarily reflects the helium defined temperature, thus protecting 
the abundance from any significant bias. 

Thus, the new MCMC method is a distinct improvement, resulting in a statistically more 
accurate determination of the helium abundance, the physical parameters associated with the 
HII region, and their uncertainties. Nevertheless, we found relatively large uncertainties in the 
helium abundance determinations of individual low metallicity HII regions. This, however, 
is an indication of the true uncertainty in the measurement and the challenge posed. Future 
work will investigate the possibilities for improving the result through the use of a revised 
and expanded set of objects. 
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